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Abstract 



We develop a model for speciation due to postzygotic incompatibility generated by autoim- 
mune reactions. The model is based on predator-prey interactions between a host plants and their 
pathogens. Such interactions are often frequency-dependent, so that pathogen attack is focused 
on the most abundant plant phenotype, while rare plant types may escape pathogen attack. Thus, 
frequency dependence can generate disruptive selection, which can give rise to speciation if dis- 
tant phenotypes become reproductively isolated. Based on recent experimental evidence from 
Arabidopsis, we assume that at the molecular level, incompatibility between strains is caused 
by epistatic interactions between two proteins in the plant immune system, the guard and the 
guardee. Within each plant strain, immune reactions occur when the guardee protein is modi- 
fied by a pathogen effector, and the guard subsequently binds to the guardee, thus precipitating 
an immune response. However, when guard and guardee proteins come from phenotypically 
distant parents, a hybrid's immune system can be triggered by erroneous interactions between 
these proteins even in the absence of pathogen attack, leading to severe autoimmune reactions 
in hybrids. Our model shows how phenotypic variation generated by frequency-dependent host- 
pathogen interactions can lead to postzygotic incompatibility between extremal types, and hence 
to speciation. 



1 Introduction 



Understanding the origins of diversity is a central theme in evolutionary biology. In general, evo- 
utionary dive r sificat i on can be desc r ibed a s temporal modification of phenotype distributions 



dCovne & Orj (120041) . 



Doebeli et al. 



(|2007h ). A single species typically corresponds to a uni- 



modal phenotype distributio n, whereas this distri bution becomes bimodal (or multimodal) once 



diversification has occurred dPoebeli et al. 



(l2007h ). Thus, speciation can be described as a split- 



ting of the ancestral unimodal phenotype distribution into a two or more descendant peaks, with 
each peak corresponding to an emerging species. Traditionally, diversification and speciation 
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are thought to occur when different phenotypes are selectively favoured in different and isolated 
geographical regions, so that an initially uniform ancestral population that is geographically dis- 
persed would develop different phenotypic modes corresponding to the phenotypes favoured in 
different locations. More recently, processes of adaptive speciation, which unfold in the absence 
of geographical isolation and during which phenotype distributions become multimodal due to 
ecological int eractions such as compet ition for resources of predation, have received consider- 



able attention (IDieckmann et al. 



(120041) ) ■ In sexual populations, adaptive speciation requires that 
the different clusters in the phenotype distribution, corresponding to the newly emerging species, 
are separated by barriers to gene flow. That is, to prevent mixing between the emerging species, 
successful reproduction should occur predominantly within each emerging phenotypic cluster, 
and not between individuals with distinctly different phenotypes. The type of reproductive iso- 
lation that has been considered in most theoreti cal models of adaptive speciat i on is isolation due 



to asso rtative mating, i.e., prezygotic isolation (|Dieckmann & Doebelil (|l999h . 



Dieckmann et al. 



(120041) '). In this paper, we investigate the potential role of postzygotic isolation for adaptive speci- 



ation by cons idering models in which postzygotic isolation is caused by autoimmune responses. 



Recently, 



Eizaguirre et al. 



(|2009|) have pointed out that the immune system may play an im- 
portant role in processes of diversification. These authors primarily considered the role of the 
vertebrate immune system for prezygotic isolation (so that mating partners would be chosen 
based on immune system characteristics), but they also mentioned that the immune system may 
be important for postzygotic isolation due to a weakened immune response in hybrids. Here we 
consider a different type of hybrid disadvantage due to malfunctioning of the immune systems, 
which is is based on observations made in plant systems. In plants, a well-known example of 
postzygotic isolation is hy brid necrosis, defined as a set of highly deleterious and often lethal 
phenotypic characteristics (^Bomblies & Weigell (|2007|) ). In hybrid necrosis, a mixture of genes 
from different strains becomes deleterious even though the contributing genes were harmless, or 
even beneficial, in the parents. Recent experimental evidence suggests that hybrid necrosis in 



plants can be caused by an epistatic interaction of loci controlling the immune response to attack 
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by pathogens (|Bomblies et al. 



(120071) ). In this form of hybrid necrosis, inviability is caused by 
inappropriate activation of the plant immune system in the absence of pathogens. Among several 
mechanisms of pathogen recognition by a plant host cell, interactions between tw o different types 
of hos t proteins, "guard" and "guardee" proteins, are thought to play a key role (|Jones & Dangl 
(120061) ') ■ When a pathogen attacks a host cell, it often injects effector proteins that manipulate 
target proteins in the host cell and thereby contribute to the success of the pathogen. These 
host targets (the guardees) are "guarded" by other host proteins (the "guard") that monitor the 
guardee's molecular structure. When a pathogen effector induces changes in the molecular struc- 
ture of the guardee (creating "pathogen-induced modified-self"), these changes are rec ognized 



by the guard proteins, which then activate the immune response (IJones & Dangll (120061) '). Thus, 
on the one hand successful pathogen attack requires fine-tuning of the effector to the guardee, so 
that a pathogen with a given effector repertoire can only successfully attack a certain range of 
host cells (i.e., those with the "right" types of guardees). On the other hand, efficient immune 
response requires fine-tuning of the guardee and the guard, so that the guard selectively recog- 
nizes only those guardees that have been modified by a pathogen. However, if mating between 
different host strains leads to hybrids in which guard and guardees come from lineages with dif- 
ferent evolutionary paths, the guardee might recognize the guard as modified even in the absence 
of pathogen attack, which could lead to immune response and subsequent necrosis of the hybrid 
even in the absence of any pathogen. Experimental evidence of two-locus epistatic interactions 
for hybrid necrosis, of the autoimmune nature of the deleterious phenotype, and of an increased 
disease resistance of rescued from necrosis hybrids all support the hypothesis that guard-guardee 
interactions are involve d in hybrid necrosis in the well-studied model plant Arabidopsis thaliana 
(IBomblies et all (|2007|) ). A large number of hybrid necrosis cases sharing phenotypic similar- 
ities with the Arabidopsis cases indicate that this may be a common mechanism operating in a 
wide range of plant species. 

Due to the potentially high selection pressures exerted by pathogen attack, genes controlling 
immune responses are generally thought to be fast evolving. This can in turn generate strong se- 
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lection pressures on pathogens, which can lead to coevolution. In addition, such host-pathogen, 
or more generally, predator-prey interactions, are often frequency-depende nt, and it is known tha t 



this frequency-dependence ca n generate disruptiye selection in the host (IDercole et al. 



(120031) . 



Doebeli & DieckmannI (llOOOt) ). Here we present a mathematical model of adaptive speciation 
in a host plant in which diversification is by driven host-parasite interactions, and postzygotic 
reproductive isolation is caused by hybrid necrosis. The barriers to gene flow between emerg- 
ing phenotypes are generated by detrimental autoimmune reactions in individuals containing 
pathogen-resistance genes from different clusters. Specifically, reproductive postzygotic iso- 
lation is due to genetic incompatibility between guard and guardee proteins of phenotypically 
distant strains. In the following we develop a simple representation of the evolution of the guard 
and guardee proteins driven by selection for escaping recognition of the guardee by pathogen 
effectors. As we will show, the frequency-dependent selection imposed by the pathogen leads to 
the emergence of distinct strains that are reproductively isolated due to genetic incompatibility in 
the immune response. The emerging host strains correspond to two distinct paths of coordinated 
evolution of guard and guardee proteins. Mating between these strains can result either in hy- 
brids that are very disease-prone due to a defective immune system, or in individuals that show 
hybrid necrosis due to autoimmune reactions. 



2 Model description 



2.1 Phenotype space 



Several experimental observations from 



Bomblies et al. 



(l2007h are essential for the definition of 



our model. First, it has been uncovered that the autoimmune reaction plays an essential role in 
lethality of the hybrid phenotype, and that the surviving hybrids exhibit increased pathogen re- 
sistance. Second, it was shown that epistatic interactions between two loci are both a necessary 
and a sufficient condition for hybrid necrosis. And finally, it was observed that an increase in a 
habitat temperature from 16° to 23° rescues the hybrid. All these phenomena strongly indicate 
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the possibility that hybrid necrosis is caused by erroneous binding between guard and guardee 
proteins coming from different parents. In normal plants, such binding only occurs if the guarded 
protein is modified by the pathogen effector. However, in a hybrid such binding could happen 
in the absence of pathogen effectors since the guard and guardee proteins evolved independently 
and could thus have acquired a propensity for binding, i.e., a higher binding energy, without 
any additional alterations through pathogen effectors. An increase in the ambient temperature 
weakens any binding, thus decre asing a chance to provo ke the unwarranted immune response in 



a hybrid, exactly as observed in 



Bomblies et al. 



(120071) . An estimate showing that the increase 



in ambient temperature applied in the experiments can indeed cause a noticeable shift in binding 
equilibrium of guard and guardee proteins is presented in Appendix A. For our model, we en- 
visage biochemical binding between guard and guardee proteins to be the mechanism for both 
immune and autoimmune responses. 

We assume that the evolution of guardee and guard protein is described by two phenotypic 
coordinates, g and r. Each host plant individual is represented by a point in this two-dimensional 
phenotype space, with the (^-coordinate describing the state or genetic makeup of the individual's 
guardee protein, and the r-coordinate describing the state or genetic makeup of the individual's 
guard protein. Mutations in guard and guardee genes cause the corresponding point to shift in 
(r, g) space. The density of plant individuals with a particular form of guard and guardee proteins 
is described by a density distribution function h{r,g). An almost homogeneous population with 
little genetic variation in immune proteins is described by a density distribution with a single 
narrow peak, while a population consisting of several strains with different guard and guardee 
proteins is described by a multimodal density distribution. 

In our framework, the pathogen is characterized by a single coordinate, e, which describes 
the genetic makeup of the parasite's effector proteins. The effector proteins interacts with the 
guardee protein of a ho s t plan t, and we define e, as in more traditional predator-prey models (e.g. 



Doebeli & DieckmannI (120001) '). so that the pathogen attack is most effective for small distances 



|e — (7I between the pathogen phenotype e and the host's guardee phenotype g. Thus, it is in the 
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interest of the host to have a guardee phenotype that is different from the phenotype of the most 
common pathogen, and this is the mechanism that generates frequency-dependent selection. 

In the plant host, the relative values of the r and g coordinates determine the immune reaction 
of the plant to a parasite effector, as well as autoimmune reactions. We assume that there is a 
trade-off between mounting an efficient immune response in the presence of pathogen effectors, 
and being prone to deleterious autoimmune reactions in the absence of pathogens. Specifically, 
we assume that when g >> r, the probability that the guard binds to the guardee protein and 
triggers an immune response is small, independent of whether pathogen effectors are present or 
not. In this case, the plant immune system is less sensitive, making the plant more susceptible to 
pathogen attack, but less prone to autoimmune reactions. Conversely, when r >> g, the guard 
protein has a high propensity to bind to the guardee. In this case, the plant immune system is more 
sensitive, making the plant less susceptible to pathogen attack, but more prone to autoimmune 
reactions. This param etrization of phenotype space is in an accord with the observation made in 



Bomblies et al. 



(|2007|) that necrotic hybrids, being rescued by elevated ambient temperature, are 
very effective in suppressing pathogen attack. 

Quantitatively, we assume that the probability for an (r, (7) -plant to die from parasitic infec- 
tion once attacked by a pathogen is proportional to ex-p[(g — r)/(T/], while the probability to die 
from autoimmune reactions is proportional to exp[(r — g)/aAi], where a/ and aAi are system 
parameters. As a consequence of this trade-off, plants tend to survive best if their phenotypes 
satisfy r g case. 

2.2 Plant and parasite evolution 

To describe the co-evolution of the plant and pathogen populations, we use a predator-prey- 
style model for the dynamics of phenotype distributions in both the host plant and the pathogen, 
with specific terms that describe the plant-parasite conflict and autoimmune reactions. Math- 
ematically, the model is a system of two integro-differential equations that gives the temporal 
evolution of the host density distribution h{r, g) and the pathogen density distribution p{e). In 
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the integrals, the quantity h{r, g)drdg is the density of plants with guard and guardee phenotypes 
in the coordinate interval (r, r + dr) and {g, g + dg), and the quantity p{e)de is the density of 
pathogens with effector phenotype in the coordinate interval (e, e + rfe). The equations deter- 
mining the rate of change in the distributions h{r, g) and p(e) consist of birth and death terms, 
describing the increase and decrease in plant and parasite population densities due to the various 
biological components and interactions, as follows. 

• For a given plant phenotype {r^g), the rate of attack from the pathogen, and hence the prob- 
ability of infection, is proportional to a weighted sum over all pathogen phenotypes, with 
the weights reflecting how well a pathogen phenotype e can attack a plant with guardee 
phenotype g. The weight function, or "attack kernel", is assumed to be of Gaussian form 
and given by 



Ga{e - 5) = exp 



2^2 



(1) 



reflecting the fact that infection is easiest if the pathogen phenotype is very similar to the 
guardee phenotype (on some appropriate scale). Accordingly, for the given plant pheno- 
type (r, g), the probability of being attacked is proportional to 

p{e)Ga{e- g)de, (2) 

where p{e) is the pathogen density distribution. Once attacked, the probability that the 
plant individual dies is determined by the immune response (as described above) and is pro- 
portional to exp[(g — r)/(j/], reflecting how well the guard protein recognizes the changes 
in the guardee protein induced by the pathogen attack. Overall, this leads to a total death 
rate of plants with phenotypes (r, g) due to pathogen attack given by 

-5ih{r, g) exp ( - — - ) / p{e)Ga{e - g)de. (3) 



Here 5i is a constant of proportionality, and h{r, g) is the density of host plants with phe- 
notype (r, g). The minus sign indicates that this is a death term. 
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In the plant population, death also occurs due to autoimmune reactions, whose magnitude 
in plants of phenotype (r, g) is proportional to exp[(r — g)/ cfai\- The resulting death rate 
is given by 

-5,Mr,5)expf^y (4) 



For simplicity, we assume that the rate coefficient 5i for the autoimmune and parasite- 
induced death terms is the same. 

To complete the death terms for the plant, we assume that density dependent competition 
in the absence of the pathogen results in a logistic death term of the form 

h{T,g)H{t) 
K{r, g) 

here H{t) is the total density of the plant population, i.e. 

H{t)= j h{r,g)drdg. (6) 

The parameter 5c in eq. (5) is a rate coefficient for the plant death rate due to competition. 
The function K{r, g) is the carrying capacity function, which we assume to be of the form 

{r-rof^ \ (g-goV 



K{r, g) = exp 



exp 



(7) 



2(72 

This reflects the assumption that due to costs and benefits of the guard and guardee that 
are unrelated to their effect on immune response, there are optimal values tq and go for 
these traits. This means that in the absence of pathogen attack and immune response, the 
plant traits would evolve to their optimal values ro and go (and the plant distribution would 
converge to a narrow unimodal distribution centred at (ro, go))- In this way, the carrying 
capacity function K provides a component of stabilizing selection that is independent of 
host-pathogen interactions. 

In general, it is not clear what the position of (ro, go) would be in phenotype space. In 
particular, it is not clear what cost stabilizing selection would impose on a well-tuned 
immune response, for which r ^ g. Therefore, we assume in the following that only the 



trait associated with the guardee (g) protein affects the carrying capacity. As this is also 



the trait that affe cts the host-pathogen interact ion, this corresponds to commonly made 
assumptions (e.g. Doebeli & Dieckmann (2000.) ). Mathematically, that fact the the r-trait 
does not affect the carrying capacity corresponds to the assumption that the width of the 
carrying capacity in the r-direction is very large, cTf. — oo. 

To derive the birth term for the plant population, we need to incorporate sexual reproduc- 
tion, and for simplicity we assume that individuals are haploid and have two loci with 
continuously varying alleles encoding the two phenotypes r and g. An offspring with phe- 
notype (r', g') either inherits the two alleles r' and g' from different parents, or from the 
same parent. In the first case, the offspring comes from a mating between (r', g") and 
(r", g') (or vice versa) with all possible g" and r". Such matings occur with probability 
proportional to '%H{t) ' where H{t) is the total plant density as before. In the second 
case, the probability that such an offspring is produced is simply proportional to h(r', g'). 
Adding the two cases together, the total probability that an (r', g') offspring is the result of 
mating is thus 



h{r',g")h{r",g'l^^„^^„^h{r',g') 



(8) 



2H(t) ' '-^ ' 2 
(the factor 1/2 reflects ambiguity in assigning mother and father in any given mating pair). 
In addition, we assume that mutation in the plant traits is described by two normal mutation 
kernels 

(r — r')^ 



1 



and 



GM,gi9 - a') 



271(7 M.. 



1 



exp 



(9) 



exp 



l\2 



-9') 



(10) 



'2710 M,g 

Here GM,r{r — r') describes the probability that a mutation in an r'-allele produces a value 
in the interval [r, r + dr], and similarly for GM,gig — g')- Thus, small mutations are more 
likely than large ones. 
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Overall, the rate at which offspring with phenotype (r, g) are produced is then given by 




P / / GmA^ - r')GM,g{9 - 9') 




h{r',g")h{r",g') ^ „ h{r',g') 
2H{t) ^5 + 2 



dr'dg\ 
(11) 



where (5 is the birth rate. 



• For the pathogen, birth is determined by successful infection. For pathogen type e', the 
probability that it successfully attacks and subsequently infects a host of type (r, g) is 
proportional to 



/i(r, g) exp ( ^^^^ Ga{g - e'), 



(12) 



where Ga is the attack kernel given by eq. (1). Therefore, the total probability for pathogen 
type e' to successfully infect any host plant is 




9 — ^ 

/i(r, g) exp ( ) Ga{g - e')drdg. 



(13) 



For simplicity, we assume that reproduction is asexual in the pathogen, but as in the plant 
species we include mutation given by a normal function GM,e (with a width described by 
a parameter aM,e, so that the total rate of production of offspring of type e is 



aSi / GM,e(e - e')p(e') 




h{r, g) exp ( ^—^ ) Ga{g - e')drdg 



de'. (14) 



Here the conversion coefficient a indicates how many new parasites are produced from an 
infected host. 

• Finally, as in many other predator-prey models the death rate of the pathogen is assumed 
to be 



-8pp{e), 



(15) 



for some parameter 5p describing the intrinsic per capita death rate of the pathogen. Note 
that we assume that this death rate is independent of the pathogen phenotype e. 
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Collecting all the birth and death terms into two equations describing the dynamics of the plant 
and pathogen density distributions now yields the following system of coupled partial differential 
equations: 

h{r',g")h{r",g')^ „^ „ , h{r',g') 



dh(r,g) 



dt 



- (^/^(^, g) exp {^-^^ j P{^)Ga{e - g)de (16) 

- 5ih{r, g) exp - 8c— rr, 

\ ai J K{r,g) 



dp{e) 



ol5i j GmA*^ - e')p(e') J J h{r,g) exp ^ ^^id - e')drdg 



de' 



dt 

-Spp{e). (17) 

These nonlinear equations have many parameters and in principle may exhibit a variety of 
dynamic regimes. Due to the apparent complexity of the dynamical system, we do not expect any 
analytical results to be feasible and instead investigated this system using numerical simulations. 
We are particularly interested in those regimes that lead to multimodal equilibrium distributions 
in the host plant. 

Such dynamics are illustrated in Figure 1 and occur for a substantial range of parameters and 
correspond to cases in which the host-pathogen interaction leads to plant speciation. Intuitively, 
the choice of parameters generating this scenario can be explained as follows: The width of the 
carrying capacity ag defines the phenotypic space scale of the model and is taken to be one. 
One of the birth coefficients, defines the time scale and is also taken to be one. Finally, the 
competition death term defines the third scale, the amplitude of the population density and is 
taken to be one as well. The mutation width for the host crM, r, aM,g and pathogen aM,e define 
the minimal width of the emerging pattern and should be significantly smaller than one. To enable 
pattern formation due to host-parasite interaction, the parasite attack width should also be less 
than one. Finally, the intensity of the autoimmune reaction and the susceptibility to the parasite 
infection should grow sufficiently fast (faster than carrying capacity decay to be relevant) as the 
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Figure 1: Left panel: Contour plot of the host population density in logarithmic scale. Clearly visible 
along the r = g diagonal are two peaks, centered around two coordinates {ri,gi) and {r2,g2), 
corresponding to two strains, driven into separation by parasite-plant interactions. Two much 
weaker peaks with coordinates approximately {ri,g2) and (r2,gi), correspond to hybrid het- 
erozygotes in the R and G genes. The subpopulation of necrotic hybrids lies below the r = g 
diagonal, where r > g. The plot corresponds to the solution of eqs. (16) and (17) at a steady 
state (t = 800) for the following parameter values: P = 1, a = 1, di = 5, dc = I, Sp = 1, 
(7/ = 0.2, aAi = 0.2, fja = 0.8, aM,r = (^M,g = CFM,e = 0.1, and ar = CFg = 1. To avoid 
unrealistically high autoimmune and pathogen-induced host death rates, they were truncated by 
replacing expit (^-^) by min [expit {J-^) , 10^]. Right panel: Plot of the host population 
/i(r, r) along the r = 5 diagonal (black solid line) and the parasite population p{e) (red dashed 
line). 
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phenotypic distance from the optimum, r ^ g, increases. Thus, the immune and autoimmune 
reaction widths aj and cr^/ must be noticeably less than one as well. Other rates are set equal to 
unity except for the coefficients for the autoimmune and parasite-induced death term, 5i, which, 
to make this term more significant, is set equal to 5. 

In Figure 1, each of the main maxima in the host plant distribution, located at {ri,gi) and 
{^2, ) 5 respectively, corresponds to an emerging host strain or sub-species with a distinct genetic 
makeup of guard and guardee proteins. Note that ri gi and r2 ~ g2, which means that in each 
of the corresponding strains, the function of the guard and guardee proteins are geared towards 
both efficient immune response to the pathogen and low likelihood of autoimmune reactions. 
The equilibrium distribution in Figure 1 also shows two secondary peaks located approximately 
at (ri, (72) and (r2, gi). These peaks correspond to hybrids between the two emerging subspecies. 
These hybrids inherit their guard and guardee genes from parents of different subspecies and thus 
are either easy victims of pathogen attack (when g > r) or exhibit strong autoimmune necrosis 
(when r > g). Despite the equal per capita rate of birth of the homozygote and heterozygote 
offspring (see eq. (11)), the peaks corresponding to the hybrids are much smaller due to their 
much higher death rates. This is illustrated in Figure lb, which shows the bimodal equilibrium 
density distributions of the host (continuous line) and pathogen (dashed line) along the diagonal, 
as well as the trimodal density distribution of the host along the anti-diagonal (stippled line). Here 
the mode in the middle corresponds to the saddle between the two modes of the host distribution 
along the diagonal. 

Obviously, the model can also produce other equilibrium distributions, depending on the 
values of model parameters, particularly the parasite attack width. Examples are given in Figures 



2 and 3. In Figu re 2, the distributions converge to a unimoda 



previous results (|Dercole et all (|2003|) . 



equilib rium. In accordance with 



Doebeli & DieckmannI (|2000|) ). this tends to happen for 



larger predator attack widths 0"^. In contrast, for small attack widths, the equilibrium distributions 
may have more than two modes along the diagonal, and hence more than two secondary hybrid 
peaks, as illustrated in Figure 3. 
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Figure 2: Same as Fig. 1, but for a larger parasite attack width, ctq = 1.2. In this case the host density 
does not evolve to a multimodal distribution, but forms a broad single peak, and the same is true 
for the pathogen density distribution (not shown). 




Figure 3: Same as Fig. 1, but for a smaller parasite attack width aa = 0.4. In this case the host density 
evolves to a multimodal distribution, with multiple secondary peaks corresponding to hybrids 
(Fig. 2a). The multiple peaks are also evident if the host distribution is restricted to the diagonal, 
and the pathogen distribution also has multiple peaks (Fig. 2b). Along the anti-diagonal, the host 
density distribution reflects the multiple secondary hybrid peaks (Fig. 2b). 
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3 Discussion 



Both 



Bomblies & Weigell (l2007h and 



Eizaguirre et al. 



(120091) recently argued that the immune 



Eizaguirre et al. 



mm 



system could play an important role in speciation processes. While 
mainly considered the vertebrate immune system as a potential source of prez ygotic reproductive 



Bomblies & Weigell (|2007h argued that defec- 



isolation after divergent adaptation to pathogens, 
tive autoimmune reactions in hybrids could be the source of postzygotic isolation. They were 
able to support t his pe rspective by impressive experimental evidence in Arabidopsis thaliana 



(IBomblies et al. 



(|2007|) '). Our aim here was to provide mathematical evidence for the feasibility 
of adaptive speciation due to host-pathogen interactions when isolation is caused by postzygotic 
autoimmune deficiencies. 



Bomblies et al. 



(|2007|) have shown that hybrid necrosis, defined as a set of deleterious and 
often lethal phenotypic characteristics, can be caused by an epistatic interaction of loci con- 
trolling the immune response to attack by pathogens in Arabidopsis thaliana. More precisely, 
hybrid necrosis is caused by detrimental activation of the plant immune system in the absence 



of pathogens, and 



Bomblies et al. 



(|2007|) showed that epistatic interactions between two loci 



are both necessary and sufficient for this form of hyb rid necrosis. This corresponds to classic 



Dobzhansky-MuUer incompatibilities (iGavrilets 



(120041) ). which we have modeled here by assum- 



ing that the optimal compromise between the ability to respond to pathogen attack and the avoid- 
ance of undesired autoimmune reactions is attained by individuals with genotypes {r,g ^ r), 
where r and g represent two loci controlling the immune response. Thus, two genotypes (r, r) 
and (i?, R) with R very different from r can both be optimal, but their heterozygous hybrids 
(r, R) and (i?, r) will suffer from a malfunctioning immune system. In the case of the plant 
immune system, the two genes r and g represent two different proteins, "guard" and "guardee" 



proteins. Such proteins are thought to be central for the plant immune system (|Jones & Dangl 
(|2006|) ). and their optimal functioning requires concerted evolution at both loci. 

The model presented here combines a detailed, albeit schematic, description of the plant 
immune response based on the traits r and g, with the macroscopic, population-based represen- 
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tation of ecological interactions driving the evolution of these traits. The basic result is that for 
a range of intuitively appealing parameters, the model leads to adaptive speciation: as a result 
of frequency-dependent selection exerted by pathogen attack, a single ancestral host strain splits 
into two descending strains (r, r) and {R, R) that each have their distinct genetic makeup of the 
guard-guardee pathogen recognition system. Hybrids with guard and guardee genes coming from 
different parents either suffer necrosis due to autoimmune reaction, or they exhibit a weakened 
immune system due to a compromised ability to detect the parasite attack. This latter form of 
hybrid disadvantage corresponds to the m echanism of postzygotic isolation due to MHC -based 
immune responses that was conjectured by lEizaguirre et ah (t2009^) to operate in vertebrates. 



In contrast to most previous models of adaptive speciation, in which reproductive isolation is 
based on prezygotic mechanisms such as assortative mating, in our models reproductive isolation 
between diverging lineages emerges due to postzygotic hybrid inviability. Of course, a natural 
question would be whether such postzygotic isolation would select for prezy gotic isolation in 



the fo rm of assortative mating based on the immune system, as envisaged by 



Eizaguirre et al. 



(|2009|) . Our models could be extended to pursue this, as well as a number of other extensions. 
For example, it is known that many pathogens are capable of producing a number of different 
effector proteins that allow them to attack a host plant. Similarly, the immune system of the plant 
has a number of different ways of dealing with this effector variety. Accordingly, it would be 
interesting to see whether the type of diversification observed in our models would be easier or 
harder to obtain if the dimensionality of both the guard and guardee traits in the host and the 
effector trait in the pathogen is increased. For example, it is possible that hybrid necrosis, and 
hence postzygotic isolation, is increased if more than one guard-guardee pairs are involved in 
the immune response to a particular pathogen, because simultaneous incompatibility of various 
guards and guardees could lead to more severe autoimmune reactions. 

Instead of considering larger numbers of guard-guardee pairs, another way to increase the 
the complexity of the model would to develop a more detailed, mechanistic description of the 
genetic network of activation and repression of the various pathways involved in the immune re- 
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sponse based on a single guard-guardee pair. In the present model, this network is assumed to be 
extremely simple in that it is assumed that the binding affin ity of guard to guardee depends on the 



genetic distance {r — g). 



Ten Tusscher & Hogewegl (|2009|) have studied more traditional models 



for adaptive speciation based on resource competition under the assumption that the regulatory 
network determining the phenotypes important for competition are much more complicated, and 
realistic, than commonly assumed in such models. One of the main conclusion of these authors 
is that genetic complexity facilitates diversification and speciation, and it would be interesting to 
see whether similar conclusions would be reached if more genetic complexity would be incorpo- 
rated into the models presented here. 

Despite being schematic and minimalistic, our model correctly reflects some of the main ex- 
perimental observation s of postzygotic isolation due to hybrid necrosis in Arabidopsis thaliana 



(IBomblies et al. 



(|2007[) . see also Appendix). It also predicts a class of hybrids with weakened im- 



mune response to path ogen attack, whose exist ence should be investigated in future experiments. 



Overall, we agree with 



Eizaguirre et al 



(|2009|) that studying the role of the immune system for 
speciation processes is very promising. Adaptation in the immune system as a co-evolutionary 
response to pathogens may be a potent mechanism for diversification based on both prezygotic 



and postzygotic reproductive isolation. 



4 Appendix 

Here we provide an argument for the fact that an increase in ambient temperature from 16" to 23° 
can indeed cause a biologically noticeable reduction in autoi mmune reaction and he nce rescue 



Bomblies et al. 



(2007). We make 



otherwise necrotic hybrids, as observed in the experiments of [ 
a simple estimate of how a change in temperature affects the Law of Mass Action governing 
the binding-unbinding equilibrium between the guard and guardee proteins. For simplicity we 
assume that both protein can exist only in two forms, free and bound to each other forming a 
dimer. we denote the free forms by R and G, and the dimer by RG. In the limit of weak binding 
(large dissociation constant k), i.e., when both proteins are mostly in free form, the concentration 
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of the dimer is 



[RG] = [RUG]o/k, 



(Al) 



where [R]o and [G]q are the total concentrations of proteins R and G. The temperature depen 
dence of the dissociation constant is usually given by the Arrhenius form, 

E " 



k = ko exp 



KT 



(A2) 



As a consequence, the relative decrease in the concentration of RG caused by the increase in 
temperature by AT is 



[RG] 



T+AT 



exp 



EAT\ 



(A3) 



[RG]t "V KT^ 

A cell can typically recognize a change in conc entration larger than 2 0% (smaller shifts in con 



centrations are apparently perceived as "noise". 



Newman et al. 



(|2006|) ). Thus, for a very reason- 



able value of protein-protein dissociation energy E ^ 8kT, an increase in the temperature from 
16° to 23° can produce a biologically meaningful decrease in concentration of the RG dimer. It 
is at least plausible that such a decrease can eliminate improper binding between R and G, and 
hence unwanted autoimmune reactions. 
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